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I. INTRODUCTION 


A. GENERAL 

Since the early 1960s the U.S Navy has been interested in 
Unmanned Underwater Vehicles(UUV), which include Autonomous 
Underwater Vehicles (AUV) and Remotely Operated Vehicles (ROV). 
The latter are tethered to a supporting platform, where an 
operator provides necessary control signals. In this case most 
of the mission planning and on-line decision making is 
accomplished by the operator. Increasing computer capabilities 
have decreased the need to tether the vehicle to prove 
guidance and navigational support. The improvements allowed 
underwater vehicles to be autonomous. [Ref. 1] 

Autonomous vehicles can go where humans cannot go and 
do not want to go. Autonomous vehicles are capable of 
receiving initial input, moving to another location, executing 
a mission, and returning with the requested results and the 
data. In addition to performing labor intensive and repetitive 
tasks, these vehicles can perform their tasks faster and with 
greater precision than humans, and can also proceed into 
hostile or contaminated environments. [Ref. 2] 

At the Naval Postgraduate School the interest in 


autonomous underwater vehicles has manifested itself in 


current research which supports development of a prototype 


vehicle named NPS AUVII. 


B. GOAL 

If an autonomous underwater vehicle is to be completely 
autonomous, it must have a navigation system capable of 
estimating its own position at all times. The effort in this 
thesis is mainly concerned with the navigation problem of the 
AUV. Compared to autonomous land vehicles, AUVs have many 
disadvantages, such as the presence of highly nonlinear 
vehicle dynamics that are often uncertain and the presence of 
currents whose effects cannot be detected by inertial 
navigation systems. One possible solution to this problem is 
the use of velocity measurements relative to the ocean floor, 
using, for example, a doppler sonar. However, the cost, size 
and the depth of the ocean often prevents use of doppler 
sonar on standard AUV’s. Visual information (sonar and 
optical) seems a viable approach in environments having 
landmarks as references, such as buoys, piers, rocks, and man 


made structures. 


C. METHOD OF APPROACH 

This thesis develops a short range navigation algorithm 
for the planar motion of the NPS AUVII vehicle. The system 
uses sonar returns and incorporates them with a potential 


function which defines the map of the operation area. The 


measurements obtained from sonar are filtered through a Kalman 
estimator uSing extended Kalman filtering. The main feature of 
the algorithm is the use of a potential function to define the 
area of operation. By use of this potential function, the 
System discriminates the objects which are not on the map , 
and yields a robust and reliable solution to the short range 
navigation problem of the AUV. The algorithm is first applied 
in a rectangular shaped pool and later applied to an 
environment without any borders (e.g open sea). 

This thesis consists of five parts. Chapter II discusses 
the theory behind the navigator design and the modeling 
assumptions in this study. Chapter III discusses the details 
of the potential function approach in the implementation of a 
navigation algorithm together with simulation results. The 
results using the actual data collected in the water tank 
are also discussed. Chapter IV discusses the implementation of 
neural networks for the navigation of AUV. Consequently 
Chapter V summarizes the results and conclusions of the study, 


and provides recommendations for further investigation. 


phat BACKGROUND 


A. KALMAN FILTERING 

Since it was derived by R.E. Kalman in 1960, the Kalman 
filter has been used extensively in the design of optimal 
estimation algorithms. Its rapid acceptance and subsequent 
success are due to the Kalman Filter’s recurSive nature and 
optimality. 

The Kalman filter estimates the states of a system given 
a set of inputs to the system and a set of measurements both 
corrupted by noise [Ref. 3]. The system is assumed to be 


driven by both a known input and a random input: 


X(K+1) =®x(k) +A,u(k) +A,w(k), (2.1) 


where u(k) is the deterministic input, w(k) is plant driving 
noise, @is the state transition matrix and x(k) is the state 
vector. We also assume w(k) iS a zero mean gaussian random 


vector with covariance: 


E(w(k)w(k+n)]=06(n), (2.2) 


where Q is the covariance matrix of plant noise and 6(n) is 


impulse function. The measurement process can be modeled as: 


y (k) =Hx(k) +v(k), Een 


where v(k) is zero mean additive gaussian noise and H is the 


measurement matrix. Measurement noise v(k) has covariance: 


E[v(k) v(kt+n)j]=R6(n), (2.4) 


where R is the meaSurement error covariance matrix. 
Using the orthogonality principle it can be shown that the 


recursive solution to this problem has the following form: 
Peete ea et eG ly ery ( Kae) (2.5) 


Given the current estimate we can predict ahead by simply 


using the state equation: 


R(k+1(k) =OR(k[k) +A u(k) , (270) 
which yields the optimal predictor. We can reasonably estimate 
¥(k+1) by simply using the H matrix such that: 

Yuri) = oe (ed ik) : (2.7) 


The Kalman gain G depends on the covariance matrix of 


eseimation error which is defined as: 


Jee [Ne Mc rel a bree ee [3h a G.c 6) eal) alle (2.8) 


The Kalman gains G(k+1) are found by applying the 
orthogonality principle as shown in [Ref. 3], which leads to 


the following recursive equations: 


G(k+1) = P(k+1|k)H7(HPH? + R)~ (2.9) 
P(K+1 Kt 12 Gl ee pee ni) (2.10) 
P(k+1|k) = O2(k|k)o7 + A,OA,7 (Zan) 


The filter equations can be summed up as follows: 


&(k+1|k) = O&(k|k) + A,u(k) (2.12) 
9(k+1|k) = H&(k+1|k) (2.13) 


R(K+1|k+1) =A(KtL | K) +e (ky Gey ea (2.14) 


As can be seen from the above equations, for a Linear Time 
Invariant (LTI) system, the Kalman gains can be computed off- 
line and stored in the computer and used as look-up tables. In 
order to start the filter algorithm, a priori values of the 
State estimates and the error covariance are needed. The 


proper choice for these estimates are: 


R(O|0) = E[x(0)], (2. 5) 


PO |\QieGov [=.(ciale (2.16) 


The Q matrix is the covariance of the state excitation 
noise. If the disturbances to the system are uncorrelated, 
the Q matrix is diagonal. The Q matrix accounts for the noise 
processes in the system as well as modeling errors between 
model and system dynamics. 

R matrix is the covariance matrix of measurement error. 
Large values of R means that the measurements are not 
consistent with the dynamic model. 

The matrices Q and R represent a trade-off between the 
model and the observation accuracies. In particular if the 
entries of R are larger than the entries of Q, it means that 
the model is more reliable than the observations, and the 
Kalman Filter will give less emphasis to the observation. Vice 
versa, if the entries of R are small compared to Q , it means 
that the observations are more reliable than the model, and 
the Kalman filter puts more emphasis on the observations. 

While the matrix R is in general determined by sensing 
device, the determination of Q is more difficult due to its 
lack of physical meaning, and it is usually set by trial and 
error. 

As mentioned above for the implementation of the filter 
algorithm we need a priori information for both the states and 
the error covariance. Large values in the initial estimate of 
P means that the filter will not solely depend on the initial 


conditions and gives more emphasis to the measurements. 


B. EXTENDED KALMAN FILTERING 

The Kalman filter yields a statistically optimal estimate 
for the states of a LTI system. Unfortunately many real 
systems such as the model studied in this thesis are nonlinear 
in nature. A general nonlinear system driven by both a 


deterministic and random input is: 
K(Kt1) =f (xi es) awe ae (2.17) 


where w(k) is a random forcing function with covariance Q, 
u(k) is a deterministic input and f(.) is a nonlinear 


function. The measurement process can be modeled as: 
VK) =O), Via, (2.18) 


where v(k) 1S measurement noise with covariance R and g(.) is 
a nonlinear function. 

Our task is to find the linear estimate of x(k), which 
minimizes the mean square estimation error. Solution to this 
problem 1s given by the conditional expectation. But 
generating and maximizing the conditional expectation is a 
formidable task which is-~ not suitable to on-line 
implementations. 

A simpler and of course mostly used approach is to uSe an 
extension of the Kalman filter by using Taylor series 
expansion. The resulting filter is suboptimal and does not 
give a guaranteed convergence as the Kalman filter does. In 
Spite of this,the Extended Kalman Filter (EKF) is used ina 


wide range of applications, especially in target tracking 


systems where the nonlinearity stems form the geometry 
involved in the measurement process. 
The extended Kalman filter equations are similar to the 


Kalman equations: 


R(k+1|k) =£(2(K|K) ,u(k),0,4), (2.19) 
Viktlikpe=ol eke iik) ,0,;K) ; (2210)) 
R(k+1|k+1) =R(k+1[k) +G(k+1) [y(k+1) -P(k+1[k) ] ( 229215) 


As can be seen from these equations the states are 
predicted ahead using the estimated current state and the 
known input. Since the noise is unknown, it is set to its 
expected value of zero. The Kalman gains are computed by using 
the first order linear approximation of the nonlinear model 
using Taylor series expansion. These approximations are as 


follows: 


u 
w(k) =0 


_k) ; 
Ox (k) [x (k) =8(4]i0 (2e22) 


WOOP (eC) tie) aw ek 
7 sec Lt SE A coos Saal AOE Al Shs A a I 
01k [x (49) = il I0 (2.23) 
w(k) =0 
_ Og (x(k) V(K),K) | _— 
- Ox (k) i a ae 


ee (2525) 
v(k) =0 


These linearized values are used in the gain equations: 


P. x+1|k) =@P(k|k) O7+AQA7 (2.26) 
GUK+l) PU kik) A er Rp alee (2.205 
P(k+1|k+1) =[I-G(kt1) H] P(k+1|k) (2.28) 


It is easy to see, from these equations, that the Kalman 
gains should be computed on-line since the linearized model 


depends on the current state estimate. 


C. SMOOTHING ALGORITHM 

Smoothing iS a procedure that uses all the state 
estimates and associated error covariances produced by the 
Kalman Filter and attempts to improve the accuracy of these 
estimates. 

Let k be within the time interval 0 to N, that is, OsksN, 
for some fixed N. The Kalman filter estimate at time k denoted 
by &(k+1|k) is based only on the measurements up to time k. 
The smoothed estimate is based on the measurements that 
occurred over the entire time interval. The smoothed estimate 
is denoted by &(k|N) and the smoothed error covariance is 
given by P(kIN). 

Meditch, [Ref. 4], classifies the smoothing algorithms 


into three categories: 
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Fixed interval smoothing, denoted by X(k|N) where 
k=0,1..N-1;N is a positive integer. 

Fixed point smoothing, denoted by &(k{j) where 
j=k+1,k+2,..k and k is a fixed integer. 

Fixed lag smoothing, denoted by R(k{k+L), where 
k=0,1,..N, is a positive integer, and L is the fixed time lag. 

In this study we use a fixed interval smoothing algorithm 
to smooth the estimates of the extended Kalman filter. The 


equations used in the smoothing algorithm are given below: 


A(k) =P(k+1|k+1) 67P71(k+1|k) (22.29) 
SURIN) “=-2 (et lel = ACK) (2 (eed IN) eke ky] (2.30) 


P(k|N) = P(k+1|K+1) +A(k) [P(k+1|N) -P(k+1]k)]A(k)7 (2.31) 


As can be seen from the above equations, the smoothed 
estimate is the Kalman filter estimate ®(k+1]k+1), adjusted by 
a weighted error term. This error is the difference between 
the smoothed estimate and the predicted estimate calculated by 
the Kalman Filter. The smoothed error covariance has no impact 
in the smoothing algorithm but it is a measure of how well the 
smoothing filter is working. If P(k]N)sP(k+1/k+1) this means 
the smoothed estimate is better than or equal to its filtered 
estimate. 

The fixed interval smoothing algorithm is an off-line 
procedure that operates backward in time. The smoothing filter 


1s initialized by the last filtered estimate, that is to say 


aK 


the last filtered estimate becomes the first smoothed 
estimate. It is also clear that, in order to apply the 
smoothing algorithm, the values of &(k+1/k), &(k+1/k+1), 
P(k+1/k), P(k+1/k+1) must be stored in computer prior to fixed 


interval smoothing. 


D. PARAMETER ESTIMATION 

In this section we present an algorithm to estimate the 
system parameters uSing least squares method. System 
identification is mainly concerned with fitting a dynamic 
model to the measured input and output data. Trying all 
possible dynamic models is certainly a formidable task, 
therefore we restrict our models to a certain model structure. 


We consider the following difference equation: 


A(z71) y(k) =B(z7+) u(k) +¢C(z7) e(k), (2.32) 


where e(k) is white noise and y(k) and x(k) are output and 
input sequences respectively. A, B, C are polynomials in the 


time delay operator expressed as: 


ANZ i) = tha, Zoot eee Bo 
BZ aa eee wap. (2.33) 
CUZ) 14 CZ ee GE aie 


This input output model is called ARMAX, which stands for Auto 


Regressive Moving Average Extended [Ref. 5]. A particular case 


2 


of this model is the ARX model which is used in this thesis. 
It comes from the ARMAX model with the moving average part 


removed and expressed as: 


At Zo dy Ve) SB ez ian ieee LK) (2.34) 


We can write this in regression form as: 


) =o7(k) @+e(k), (2.35) 


where 


f=[- Koi pes yp < 1 Kod eee 7 = 
Cea) Aer Rane TEENS ta oes 


The basic method to estimate the parameters @ is the least 
squares method which minimizes the mean square estimation 


eiamO : 
M 
G-argmin)’ ly(k) -$78/, (2237) 
=1 


with respect to @. The solution to this problem is given by 


pseudo inverse matrix: 


6=(©7@) 1O7y, (2.38) 


where 


ee 


(2.39) 
Sr i le 


with M being the number of data points used in the parameter 
estimation. 


If we consider a first order difference equation of the 
forme 


Y(K) +ay (hey Sai eae) K) 


(2.40) 
Then using Equations (2,38) and (2,39) 


the least squares 
estimate of the parameters a and b are given as follows: 


M 


=i 


M 
y) y(k-1)? 
K=1 








M 
=\o VR Dy ad ka sy (ko alee) 
a k=1 k=1 
2 (2.41) 
b M M M 
=) V(K=ye) Yo u(k-1)? Y> u(k-1) y(k) 
k=1 k=1 k=1 


Using either the pseudo inverse matrix form or the 


Summation form, the parameters a and b can be estimated uSing 
the above equations. 


E. SYSTEM MODELING 


Our study in this thesis is concerned with the planar 
motion of the AUV, 


therefore system modeling for the planar 


motion 1s considered. The system model used in this study is 
the same as in [Ref. 6]. 
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In modeling the planar motion of the AUV we use a earth 
fixed carteSian coordinate frame. We define the states of the 


vehicle X ¢« R as: 
oe 
C 
eee a . (2.42) 
0(t) 
Baie) 


The states x and y define the x and y position of the AUV. The 
State v is the forward velocity of the AUV and the states @ 
and § are the heading(yaw) and heading rate of the vehicle. 

The differential equations describing the model are as 
follows: 


x=vcos8 


y=vsin® 


i (2.43) 
V=-a,Vt+a,u, 


§=-b,0+b,u,, 
where a,, a,, b,, b, are constants depending on the vehicle 
dynamics and the inputs u,, u, are the RPM and rudder commands 
to the system respectively. 
The system model can now be expressed in state space form 


as 


15 








00 cos8 0 O 0 0 
00 sin6 0 0 0 0} 7, 
X =|0 0 =a OOM ee tesa Cre lea leeene (2.44) 
ceo) 0 ao 0 o| be 
00 0 0O-b 0 ge 


ray 
— 

Ny 
LH 


where W is zero mean gauSSian noise which accounts for the 
errors between the model and the actual system dynamics. 
Notice that this model does not include side-slip since we 
assume the velocity vector to be in line with the heading of 
the vehicle. For simplicity, the error due to side slip is 
taken into account by the noise term W. 

In the next chapter we are going to use this model in 
conjunction with the model of the environment to determine the 


optimal state estimator. 
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III. APPLICATION OF POTENTIAL FUNCTIONS 


A. POTENTIAL FUNCTION THEORY 

Autonomous underwater vehicles must have the ability to 
locate themselves in partially known environments using 
sensors available to them. In this study, we address the 
problem of navigating the Autonomous Underwater Vehicle (AUV) 
uSing sonar returns. These sonar returns are incorporated with 
a potential function and used as nonlinear measurements in 
Kalman Filter. The potential function is a function which, in 
a sense, defines the area where the vehicle operates. The 
basic idea behind the potential function is the following: 

1) it returns a value of zero at the boundaries of the 
object and returns a value between one and zero elsewhere; 

ii) its derivative is maximum at the boundaries of the 
object. 

As an example, if we consider a pool of rectangular shape 


with sides L,and L, we can define a function $(x,y) such that: 


Den ee ym, (3.1) 


where we take the lower left-hand corner of the pool as the 
Origin of the reference frame. If the sonar return is from the 
Sides of the pool it returns a value of zero. We combine the 
defining function with a squashing function to provide 


boundedness and define the potential function as: 


i er 


Vix, y)-= Fp ee (352) 
The function F,(z) must be smooth, bounded and differentiable. 
A possibility is to choose F,(z) to be a sigmoid functioem, 


~po-A/2 
17 (3.3) 
Lae *? 


FZ) = 
with A being a strictly positive parameter. 
The potential fFUNCEIOn = must satisfy the following 
requirements: 

1) V should be continuous and differentiable. 

11) V(x,y,.) Should return a value of zero if the return is 
from a reflecting surface. 

iii) V should be uniformly bounded. 

Note that the first and the third requirements are met by 
the sigmoid function, while the second is met by the defining 
function $8. The plot of a sigmoid function for various values 
of A is given in Figure (3.1). The function is continuous and 
bounded over the interval [-1, 1] and as the parameter A 
decreases, the "Sharpness" of the function increases. 

If we consider a field of operation containing several 


Obstacles OO) an each of which has elliptical shape, 


n? 


we can describe a defining function for each obstacle as: 


2G 
B, =[x y] A, y =6., (3.4) 
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where A ¢ R*™ and C ¢ R (Ref. 7]. The total defining function 


mo emMe product of each Lunction: 


Bix,y) = [[B;(xy). (3.5) 


: 
As an example, if we want to adda buoy of cylindrical shape 
of radius r at point(a,b) in the pool the defining function 


will be: 


B once (X,Y) = (X-a) 2+ (y-b) 2-42, (3.6) 


amgiemene total defining function will be: 


Boe (X,Y) = x(x-L,) y(y-L,) ((x-a)?+(y-b)?-4r?)), (sa) 


where the overall potential function is: 


VA ey) = ee xy) ) (3.8) 


By following the above procedure, additional obstacles may be 
added to the environment. 

The parameter A’ has a vital role in shaping the potential 
function, the small values in A result in a steeper potential 
function surface while larger values result in a smoother 
surface. In Figure (3.2) we see the 3D plot of a potential 
function in the pool. It is easily seen that at the pool 
borders it gives a value of zero, whereas at the other points 


Hem retirns a value between zero and one. 


i 


Sigmoid Function 





Figure 3.1 Sigmoid Function for various values of \. 
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Figure 3.2 Potential function for the pool. 
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B. APPLICATION TO SYSTEM MODEL 

The system model in this study has the form of Equation 
(2.44). The vehicle states are position in x and y 
coordinates, velocity, heading (yaw) and heading rate 
respectively. The navigation algorithm which is based on the 
extended Kalman Filter uses the sonar range measurements. 
Therefore we need to define a measurement equation which 
relates the states of the vehicle to the sonar range 
measurement. The measurements available to the system are 
sonar ranges p at angles aw with respect to the longitudinal 
axis of the vehicle. We define the location of the "tip" of 


the sonar range as 
[x0, yO] =(x+pcos (O+a),ytpsin (Ota) ), (3.9) 


where x, y, and @ are the current position and heading of the 
vehicle. 

We relate the sonar range measurement to the system model 
via the potential function so that the measurement equation 


becomes: 
y=g(X,p,a) =V(xt+pcos (Ota), y+psin(Ot+a) ) +W1, (3.10) 


where the function g is a nonlinear function of the states and 
the sonar measurement and W1 is zero mean gaussian measurement 
noise. The function g is the potential function which defines 


the area of operation. Apart from measurement noise the 


Za 


measurement equation should have the value of zero if the 
sonar return is consistent with the map of the area. 
By defining the above measurement equation we can write 


the overall model in the standard state space form as: 


X=AX+Bu+w 
(3.11) 


y=g(X,p,a) +W1 
The system has a nonlinear state equation and a nonlinear 
observation equation. In order to apply extended Kalman 
filtering we need to get the linearized observation matrix H 


aS in Equation(2,24), which is: 


_09(X,p,a%) _-dg dg Og 
H ae Sere a7 0, = Ole (3.12) 


The partial derivatives are found by application of the chain 


rule: 


(3°33) 


06 oz 0x 06 oz dy 0’ 


where the function F, (z) is the sigmoid function as in 


Equation(3,3). 


C. SIMULATIONS IN A RECTANGULAR POOL 
The estimation algorithm was tested first in a pool of 
Size 12X12 feet. The Tritech ST725 high resolution sonar 


onboard NPS AUVII was simulated. The sonar has a rotating head 


Ze 


with a scan rate of nine seconds and yields 400 sonar returns 
per scan corresponding to 0.8 degrees per return. The data 
from seven successive scans of the sonar head are processed by 
the Kalman filter. At point (9,9) in the pool we placed a 


buoy of radius 0.5 feet, and the potential function becomes: 


Via ee ak 2) Vy) eS) Py Oe =0225))s ate) 


fimieroeure (3.3) the position estimation with known initial 
location is shown. The location of the sonar tip, 
corresponding to the estimated borders of the pool is plotted 
for each of the seven successive sonar scans and the contour 
of the potential function is superimposed. In the last plot 
the estimated trajectory of the vehicle is shown. The 
estimated trajectory of the vehicle matches the actual 
trajectory. 

The parameter A plays an important role in the estimation 
algorithm. Its value should be kept high enough to correct for 
errors in the initial estimates and low enough to discriminate 
objects which are not on the map. The filter uses a time 
varying XA which decreases with time so that the potential 
function surface get steeper as the estimates converge to the 
actual values. When there is not enough information in the 
initial position of the vehicle the filter starts with a large 
value of A to account for the errors in the initial estimates. 

In Figure (3.4) result of the estimation with unknown 


initial position is shown. The estimated trajectory converges 


Z3 


to the actual trajectory. The effects of the time varying A 
can be seen by examining the contours of the potential 
function. Figure (3.5) shows the result of estimation where we 
have an error in the initial heading as well as the initial 
position. The system iS more sensitive to heading 
perturbations, however, the estimation converges to the actual 
trajectory. Figure (3.6) shows the result of simulation for 
the case when the obstacle at point (9,9) is not defined in 
the potential function. The uncharted obstacle is detect daa 
the algorithm and the actual trajectory is recovered from the 
estimation. 

In Figures (3.7) and (3.8) we repeat the experiment with 
two cylindrical buoys placed at locations (2,2) and (9,9) in 
the pool. In Figure (3.8) the obstacles are not defined in the 
potential function. Results of estimation show that the 
algorithm proves to be robust in the presence of uncharted 


obstacles. 
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Figure 3.3 Position estimation with known initial location, 
one obstacle present. 
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Figure 3.4 Estimation with 
obstacle present. 
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Figure 3.5 Estimation with an error in initial position and 
heading, one obstacle present. 
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Figure 3.6 Estimation with unknown initial position, the 


obstacle is not defined by potential function. 


28 









7 o 
= 2 
= = 
= > 
rm rm 
= < 
& & 
= > 
O 5 10 
X (in feet) 
7 ro 
Re ~ 
‘s & 
> > 
X (in feet) 
POOL 
ro zr 
— = 
& [s 
~ > 
0 = 10 0 5 10 
X (in feet) X (in feet) 


Figure 3.7 Estimation with unknown initial location, two 


obstacles present. 
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Figure 3.8 Two obstacles present, unknown initial location, 
the obstacles are not defined by potential function. 
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D. SIMULATIONS IN AN ENVIRONMENT WITHOUT BORDERS 

After testing the algorithm in a pool we proceed further 
to investigate its behaviour in an environment without 
borders, such as the open sea. In this case the potential 
me icon only defines the obstacles in the area of operation. 
For the case when one obstacle of cylindrical shape is 
present, the potential function will be: 


Ve BO ee) Ree) Ra a 


it returns a value of zero if the sonar range return is.from 
Beemenarted obstacles. A 3D plot of the potential function, 


for the case of one obstacle is shown in Fiqure (3.9) 


Potential Functon When One Obstacle Is present 






_cetees OY sees 
| : een is 
“rut i: 

0.6~ “AW 





0.8~ — nagar «. 
n 0.4 | ‘ 


y(in feet) wei x(in feet) 
Figure 3.9 Potential Function when one obstacle is present. 


In generating data for this simulation a high value of 
sonar range is given unless the sonar return is from an 
obstacle. We also processed the data from seven successive 
sonar scans using a time varying A parameter. Figure (3.10) 
shows the estimation when there is an obstacle of cylindrical 
shape at point (5,5) in the environment. We also assume to 
have exact information of the initial location of the vehicle. 
As we notice, the actual and estimated trajectories coincide. 
In Figure 3.11 we show the position estimation when there is 
an error in the initial position of the vehicle. In this case 
the estimated tre =ctory eventually converges to the actual 
trajectory. In the second simulation we included three 
cylindrical buoys at locations (5,5), (15,15), (18,5) "Sauemere 


the potential function description is: 


B.( x, ane (yee 0, 25 
,Y) =(x#=15)4 -15)*-0.25 
Bo(x,y) =(x ey, ) Sea 
By) =e 18) Ay Seems 


Vix, y) =F) (pate) Paegeeer ps (x, y)) . 


A 3D plot of the potential function is shown in Figure 
(3.12). The result of the estimation is shown in Figure(3.13). 
In Figures (3.14) and (3.15) we included two extra obstacles 
at points (3,10) and (8,20) which are not included in the map 
of the environment. In both cases the algorithm determines the 
obstacles that are not on the map and marks them on the 


subsequent plots. 


32 


Estimated Trajectory of AUV 
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Figure 3.10 Position estimation with known initial locations, 
one obstacle present. 
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Figure 3.11 Position estimation with unknown initial position, 
one obstacle present. 
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Potential Function When Three Obstacles are Present 
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Figure 3.13 Estimation with unknown initial position, 
obstacles present. 
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Figure 3.14 Obstacles that are not on the map are marked, 
known initial position. 
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Figure 3.15 Obstacles that are not on the map are marked, 
unknown initial position. 


a 


E. THE EFFECTS OF CURRENTS 

Currents present in the area of operation is one of the 
main contributions to the uncertainty in the position of the 
vehicle. The estimation algorithm should be modified to 
accommodate the effects of currents in the navigation 
algorithm. The inclusion of currents leads to the addition of 
two extra states in the system model. We further assume the 
Currents to be constant or slowly varying. The two extra 
aecees are C.and ©, which are associated to the components of 
the currents in the x and y directions respectively, such 


Ghat: 


(3217) 


The system states becomes X e« R’ and the system model is 


represented as: 
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In Figure (3.16) the effects of currents which cause a 
Significant amount of offset in the final position of the 
vehicle is shown. Figure (3.17) and Figure (3.18) show the 


SeGimation for known and unknown values of the initial 


Sh) 


position assuming the currents are known. We further tested 
the algorithm for the case when we do not have any information 
about the currents in the area of operation, which might 
generally be the case. Figure (3.19) gives the results of 
position estimation, which shows a reliable estimate in the 
presence of unknown currents. 

In Figures (3.20) and (3.21) we repeat the experiment for 
the case when there are two obstacles in the area. For both 
cases the navigation algorithm is able to recover the actual 


trajectory. 
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Figure 3.16 The effects of currents in the trajectory of the 
vehicle. 
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Figure 3.17 Position estimation for known initial location, 
three obstacles present. 
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Figure 3.18 Position estimation with unknown initial location, 
three obstacles present. 
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Figure 3.19 Position estimation with unknown initial position 
and current, three obstacles present. 
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Figure 3.20 Position estimation with unknown initial position 
and currents, two obstacles present. 
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Estimated Trajectory of AUV 
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Figure 3.21 Position estimation with unknown initial position, 
two obstacles present. 





For the simulations in this Chapter we have applied Least 
Squares (LS) parameter estimation techniques to estimate 
dynamic parameters a,, a), bD,, b, associated with Equation 
(2.44). 

The results of the parameter estimation are given in Table 
(3.1). The actual values are used in data generation for the 
computer simulations. The table shows the mean of the 
estimated parameters for six different simulations. The 
percentage error made in the estimates are less than one 


percent which is within acceptable limits. 
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Table 3.1 ESTIMATED DYNAMIC PARAMETERS BY LEAST SQUARES. 


eet Actual Value | Estimated percenraues 
Parameter Value — 
fa | 2.0006 0.9977 0.23 2 


© . OREO @-0010 6 


[| 000 1.0002 
2.0000 












F. STUDIES WITH ACTUAL DATA 

In this section we apply the navigation algorithm using 
actual data. The data is collected on a water tank of size 6x6 
meter. The data consists of a set of sonar returns including 
range p and bearing a. The sonar measurements are taken by a 
Tritech 725 high resolution sonar onboard the NPS AUVII. Five 
Successive sonar scans are processed by the Kalman filter. 
Fach sonar scan consists of 200 range measurememes 
corresponding to a sector of 1.8 degrees per sonar return. The 
typical sonar scan of the pool is shown in Figure (3.22). Due 
to gain problems in the sonar, the range measurements are 
excessively corrupted by noise. We see these effects in the 
sonar scan of Figure (3.22), even though there are no objects 


inside the watertank we have sonar readings inside the tank as 
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well as the outside of it. In collecting data for the 
experiment, no RPM or rudder command was given to the vehicle 
and constant speed and heading was used. Therefore, we used a 
Simpler version of the system model. In this case our state 
vector is X e« R*, given by x and y position, velocity and 


heading. This results in a state space form of 


0 

oe (ayo) 
0 

0 


where W is zero mean gauSSian white noise. In this case we do 
not have any obstacles in the tank so the potential function 


will only define the borders of the pool such that: 
VA) = Ba ee eG V6) i. (33720) 


The range and bearing values from each successive sonar 
return are fed to the Kalman Filter and the parameter A is 
reduced by 50 percent at each complete scan. The position 
estimation for known initial location iS given in Figure 
(3.23). In the subsequent plots, estimated borders of the pool 
for each successive sonar scan are plotted along with the 
contours of the potential function. Also the estimated 
trajectory of the vehicle is shown. The readings inside the 
pool are due to the gain problems associated with the sonar as 
mentioned before. In Figure (3.24) the position estimation 


with unknown initial condition is shown, where we have the 
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center of the tank, namely the point (3,73) as an dnmtecem 
estimate. The estimation converges to the actual values. 

The fixed interval smoothing, as explained in Chapter II, 
is an off-line procedure that uses the state estimates and 
associated error covariances to smooth the estimates produced 
by the Kalman Filter. In order to overcome the transient 
problems in the filter and get better estimates about the 
position of the AUV we applied the smoothing algorithm to the 
Kalman filter output. In applying the algorithm we stored the 
state estimates and the error covariances for each scan and 
fed these values to the smoothing filter. The smoothing filter 
runs backwards in time and the final state estimate and the 
error covariance of the Kalman filter become the a priori 
estimate and the error covariance of the smoothing filter. 

In Figures (3.25) and (3.26) results of the estimation 
with a smoothing filter are shown. The smoothed trajectory is 
given in the last plot. It is easily seen that fixed interval 
Smoothing provides better state estimates and does not suffer 
the transient problems associated with the Kalman filter. In 
Figure (3.27) we applied the Fixed Interval Smoothing 
algorithm to the computer generated data. The main 
disadvantage of the smoothing algorithm is having to store a 
large amount of data, which may cause memory overflow in the 


computer of the AUV studied in this thesis. 
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Sonar Scan in the Water Tank 
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Figure 3.22 Typical sonar scan in the pool. 
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Figure 3.23 Position Estimation with known initial location. 
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Figure 3.24 Position estimation with unknown initial location. 
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Figure 3.25 Fixed interval smoothing, known initial location. 
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Figure 3.26 Fixed interval smoothing, unknown initial location 
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Figure 3.27 Fixed interval smoothing applied to computer 
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IV. NEURAL NETWORK MODELS FOR POTENTIAL FUNCTIONS 


A. GENERAL 

A neural network consists of many simple elements 
operating in parallel. It can be trained to perform complex 
tasks involving system identification and classification, 
control systems, vision and pattern recognition [Ref. 8]. 

Neural networks are composed of neurons. Each neuron in 
the net is made of weighted inputs, a bias and a transfer 
function. The input to a neural network is multiplied by its 
weight and summed by the bias b. The role of the bias is to 
shift the argument of the transfer function by a corresponding 


amount. The model for a neuron is given in Figure (4.1). 


Neuron Layer 





Figure 4.1 Compact notation for the neuron model. 


oi 


In this representation R represents the number of inputs to 
the neuron and Q is the number of input vectors which is 
called batch. Two or more neurons can be combined ina single 
layer and a neural network may contain one or more such layers 
[Ref. 8]. 

When a neural network consists of multiple layers of 
neurons, the layer whose output is the overall network output 
is called output layer and other layers are called hidden 
layers. An example of a two layer neural network is given in 
Figure (4.2). In the hidden layer there are $1 neurons and in 


the output layer there are S2 neurons. 


Input Hidden Layer Output Layer 
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“7 “sixg Wy $2XQ 


F2 


| +162 
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Figure 4.2 Two layer neural network. 


a2 


The output at the hidden layer is given by: 
Al = F1(W1*P+Bl) . (4.1) 
The overall network output A is given by: 


A=F2(W2*A1+B2) 
A=F2(W2*F1(W1*P+B1)+B2), 


(4.2) 


where F2 is the transfer function of the output layer, W2 and 
B2 are the weights and biases for the output layer 


respectively. 


B. BACKPROPAGATION NEURAL NETWORKS 

The backpropagation algorithm is designed to train multi 
layered neural networks with different transfer functions to 
perform function approximation and pattern recognition [Ref. 
8]. It can be shown that a backpropagation neural network with 
at least one sigmoid layer is capable of approximating any 
Smooth function. 

The algorithm adjusts weights and biases so as to minimize 
the sum squared error between the network’s output and the 
desired signal. This procedure is done in the direction of 
the steepest descent with respect to the error which is called 
gradient descent procedure. 

Backpropagation neural networks give reasonable answers 
when they are fed with inputs they have never seen. This 


generalization property gives the convenience of training the 


po 


network with fewer input/output pairs instead of the whole 
input/output data) [Ret ele 
In the backpropagation algorithm the parameters W are 


updated as: 


~ “ OJ 
Wroy = Wie b> lwen, (4.3) 


where J is the square of the error between network’s output 
and desired signal, and yp is the learning rate. If the 
learning rate is too high it causes unstable learning, ina 
sense that the weights will diverge. If it is too small it 
requires very long training times. 

Another complication with backpropagation is the problem 
of local minima. This is caused by the fact that the error 
function is non quadratic in the weights. It is possible for 
the estimated parameters to become trapped in one of the local 
minima. This problem can be addressed by adding a momentum 
term. Momentum makes the neural network insensitive to the 
small variations in the error surface. Therefore it prevents 
the neural network’s solution to converge to a local minimum 


which causes higher steady state error. 
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C. MODELING OF POTENTIAL FUNCTIONS BY NEURAL NETWORKS 

In this thesis we use backpropagation neural networks in 
the context of function approximation. The main idea behind it 
is to get approximate values of the potential function and its 
derivatives for use in the Kalman filter. As stated in Chapter 
III, the nonlinear measurement equation is the potential 
function. The linearized H matrix is therefore made of the 
derivatives of the potential function with respect to the 


states: 


av av 


ov ov 
Ox’ oy 


H=[ , GF a6 ' 


oom (4.4) 


When the system detects a reflecting surface at point 
(x,y) in the plane, the network is trained with a "zero" 
output in accordance with the definition of the potential 
function as in Chapter III. 

The inputs to the neural network are the coordinates of 
the tip of the sonar beam which is a two element input vector. 
The potential function with its derivatives is a three element 
output vector which we use to train the network. These input 
and output pairs are nonlinearly related. We apply a 
backpropagation algorithm for our network model. The network 


architecture for this model is shown in Figure (4.3). 
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Figure 4.3 Neural Network architecture. 


As the input to the network we have the coordinates of the 
tip of the sonar beam. In layer one, nine neurons with a 
tangent-sigmoid transfer function are used. In layer two, 
Since our output size is three, only three neurons are used. 

The results of training for the potential function are 
given in the following figures. In Figure (4.4) a two layer 
plain backpropagation network is used. After 300 epochs the 
sum squared error reaches a minimum value. Figure (4.5) shows 


the results of training when backpropagation with added 
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momentum term is used. While the plain backpropagation 
converges to a local minimum, the momentum makes the network 
converge to a lower sum Squared error. 

The training time can also be decreased by using an 
adaptive learning rate by optimizing the uw based on the error. 
Figure (4.6) gives the results of training when an adaptive 
learning rate with momentum is used. The learning process is 
faster; while it takes 300 epochs to reach a minimum sum 
Squared error with momentum only, the same error level can be 
reached within 120 epochs by adding an adaptive learning rate. 
Consequently a two layer backpropagation neural network with 


adaptive learning rate and momentum is used for our application. 





Figure 4.4 Results of training with two layer backpropagation. 
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Figure 4.5 Two layer backpropagation training with momentum 
term. 
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Figure 4.6 Training with adaptive learning rate and momentum. 
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After training the neural network with input and output 
data sets as explained before, the values of the weight and 
bias matrices are stored. After this, the potential function 
is computed using the neural network. During the estimation, 
the coordinates of the sonar tip (X%,yY)), are fed to the neural 
network. The values of the potential function and its 
derivatives are taken as the network's output, and used by the 
Kalman filter. However, the training of the neural network for 
the case of unknown initial conditions of the vehicle was not 
successful due to the uncertainties in the problem. 

The results of the estimation based on the backpropagation 
neural network are given in the following figures. In Figure 
(4.7) there is one obstacle located at point (15,15) and in 
Figure (4.8) there are two obstacles located at (5,10) and 
(10,5). In both cases the estimated position follows the 
actual trajectory of the vehicle. Finally, in Figure (4.9), 
the algorithm is applied to the actual data collected in the 
water tank, which exhibits fairly good estimates about the 


D@swel10n Of the vehicle. 
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Figure 4.7 Estimation using neural networks, one obstacle 
present. 
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Figure 4.8 Estimation using neural networks, two obstacles 
present. 
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Figure 4.9 Estimation using neural networks with data 
collected in the water tank. 
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Vi SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 

In this study the problem of navigating an Autonomous 
Underwater Vehicle in a partially known environment was 
investigated. In the navigation algorithm we made use of 
potential functions to define the area of operation. The 
potential function was combined with the vehicle dynamics in 
an extended Kalman filter. Throughout the simulation studies 
the algorithm proved to be robust even in the presence of 
obstacles that are not included in the map. 

The method was first tested in a rectangular pool 
featuring several cylindrical buoys in the environment. 
Simulation studies were performed assuming known and unknown 
initial locations of the vehicle. Subsequently, the algorithm 
waS applied in an open environment, such as the open sea, 
which included several obstacles. The obstacles that were in 
the area but which were not included in the map were 
identified by the algorithm. The effects of currents were also 
sought. Simulation studies showed that even in the presence of 
currents the algorithm gave a reasonable solution to the 
problem. 

Finally, the algorithm was applied to the actual data 
collected in the water tank, The results with the experimental 
data proved that the algorithm could be reliably used in the 


implementation. A Fixed interval smoothing algorithm was also 
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used in order to smooth the estimates produced by the Kalman 
filter. Throughout the simulation studies the dynamic 
parameters were estimated by least-squares techniques and 
these estimated parameters were used in the filter. 

In Chapter IV a neural network was designed for use in the 
navigation algorithm. Basically, the potential function was 
replaced by a backpropagation neural network. This approach 
was successful only in the case of a known initial location. 

The simulation studies using both computer generated and 
the actual data show that the navigation algorithm could be 
used in NPS AUVII. Several topics might further be 
investigated such as: 

Including more complex features in the potential 
function so that a more precise map of the area could 
be used; 

e Training neural networks so that the neural network 
can be used in the presence of uncertain initial 
conditions; and 

e Converting the navigator main program to C”* 
programming language for use in the NPS AUVII on board 


GOMpUGer. 
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APPENDIX 


A. PROGRAM STRUCTURE 

The programming which covers an integral part of this 
study is accomplished using the MATLAB software package. The 
main program is NAVIGATE.M. This program reads the data file, 
prompts the user for the values of the initial state 
estimates, error covariance and the parameter A. There are 
three subroutine calls in this program. PAREST.M estimates the 
dynamic parameters associated with the system model using a 
least-squares algorithm. EST12P.M is the function that 
commands the estimator; it feeds the data to the Kalman 
filter, adjusts the parameter A and does the plotting 
routines. SCAN10.M is the function that applies extended 
Kalman filtering. VP512.M is the function that calculates the 
potential function and its derivatives in the rectangular 
pool. SIGMOID.M and DSIG.M are the functions that gives the 
value of the sigmoid function and its derivative. SHOW10.M is 
the plotting routine. EST12P.M ,SCAN10.M ,SIGMOID.M, DSIG.M 
and SHOW10.M are modified from Ref[6]. 

VP_CYL.M is the function that gives the potential function 
for the rectangular pool when there is one cylindrical 
obstacle present. EST_OPNP.M is the function that commands seme 


estimator when there are no borders. VP_OPEN.M and VP_OPEN2.M 


64 


are the functions that give the potential function in an 
environment without borders. PLOT _OPEN2.M is the plotting 
routine for the potential function. 

EST _CUR1.M and SCAN CUR1.M are functions used to model the 
estimator in presence of currents. EST _ACT2.M and SCAN ACT2.M 
are the functions that are used with actual data collected in 
the water tank. These functions also apply fixed interval 
smoothing to Kalman filter estimates. SMOOTH.M is’ the 
subroutine that accomplishes the fixed interval smoothing. 
This routine must be provided previously stored state 
estimates and predictions and associated error covariances. 

TRAIN2.M is the M file that applies backprogation with 
adaptive learning rate and momentum using Neural Networks 
Toolbox of MATLAB. The program trains the network given input 
and target output sets and generates the weights and biases 
for the neuron layers. EST _NEURO.M and SCNNEURO.M are the 
functions that are used in neural network approach. 

PR _DATA6.M and RANGEH.M are the M files modified from 
Ref[6] which are used to generate data for the simulations in 
the rectangular pool. PR _OPEN.M is the program to generate 
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Naviga 


Zim 


Estimate dynam 


e 
s 


tor algorithm 
(data, P,zh0, lambda, prms) ; 


Continue With another ysimulataen 


t(’Do you want more simulations ? Y/N,’s'); 


prms= 

=est act4 
P= eau 
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ecis (’ normal” ) 
subplot (221) 
zhm=zeros (5,1); 
ns=400; 


for i= O:n/ns-l 


fi2a., Zo) =SsCanio(zh0, P,data((ns*i+1): (ns*(i41)),:),lambda, pr 
Meal, L2) ; 
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end 
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the map and final error covariance. 

P 5 x 5 matrix,initial covarience matrix 

zho 5 x 1 matrix initial values of states 
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Ll and L2 are the sides of tthe rectangular pool. 
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R=1; 

al=prms (1) ;a2=prms (2) ;b2=prms (3) ;bl=prms (4) ; 
TS=0502257 

[Or tE— 20-1) 


rho=data(t,2); 
alpha=data(t,1); 
headh=zh(4,t); 
RPM(t) =data(t,6); 
RUDDER (t) =data(t,7) ; 
a=cos (headh*pi/180) ; 
b=sin(headh*pi/180) ; 
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[Phi,del]=c2d(A,B,Ts); 
dx0=rho*cos ( (alpha+headh) *pi/180) 
dy0=rho*sin( (alpha+headh) *pi/180) 
xO=Zh(1,t) +dax0; 

yOuzh (272) +dy0- 

Z3 (2) 6) = beOng,0 lle 
[v,dvx,dvy] =VP_cyl (x0,y0, lambda) ; 

h= [dvx,dvy,0, (-dvx*dy0+dvy*dx0) *pi/180,0]'; 
s=h’*P*h+1; K=Phi*Pth/s; 

e=0-Vv; 

zh(:,t+1)=Phi*zh(:,t)+del* [RPM(t) RUDDER(t)]‘’+K*e; 
P=Phi*P*Phi’ -K*s*K’; 
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rpm=data(:,6); 


for t=2:length (data) ; 
als=al+(v(t-1))°2; 
a2=a2+(v(t-1)*rpm(t-1)); 
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a4=a4+(rpm(t-1))*2; 
eee Wien) eit 
a6=a6 + rpm(t-1)*v(t); 


end 
% discrete velocity parameters 
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6 estimation for heading parameters 
thdot=data(:,9); 

Buaaer=Gata(:,7); 

mo 2=0 -b4=0-b5=0-b6=0- 


bom t=2:-length(data) ; 


bl=bl1+(thdot (t-1))%*2; 

b2=b2+ (thdot (t-1)*rudder(t-1)); 

b4=b4+ (rudder(t-1))%2; 

bo bS+EhGot (t-1) *thdot {t)- 

wio—O6 + rudder (t-=1) *thdot(t); 
end 
eps=le-6 
+6 discrete heading parameters. 
thhead=inv( [bl -b2;-b2 b4])*[-b5;b6]; 
phi=[-thvelo(1) 0;0 -thhead(1)]; 
del=[{thvelo(2) 0; O thhead(2)]; 
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Z=(VO0) / (lambda); v=sigmoid(z) 
Gece asig(z)*((x-12)*(y.* (y-12)) ’+x*(y.*ly- a ) /lambda; 
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end %potential 
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[v,dvx, avy] =vp cyl(x,y, lambeday) 


potential function for the pool 
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tae Lb 2=12 > 
a=9;0=9;Cc=.04; + radius of the object is .2 ft. 
vO= os * (xX- Li). Eee 2 yey ie) ee kL) eee y= ie) 
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end evp_cyl 

Beaecelon {zhm]=est opnp(data,P, zh0, lambda, prms) 

* function zhm=est_opnp (data, P,zh0, lambda, prms) 

+ modified April 1 st 1994 

6 this function commands estimator for the case of open sea 

ela ae lc: 

n=length (data) ; 

ns=200; 

moet = O:n/ns-1 

[zh,P,zs]=sen_opn(zh0,P,data((ns*i+1): (ns*(i+1)),:),lambda,p 

rms ) ; 


Znm=|zhm, zh] : 

zsm=(zsm, zs]; 
lambda=.3*lambda 

ma lengem( zh)? zh0=znis.nz) 


end 
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lambda, pause 

deasi (0 25 -5 20)) 

mieeeoon2 (25,15); 

elke. .oOn 

mies zhi l,:), zhm(2,:),"'*Gg",data(:,4),data(:,5),’-.'), 
title(’ Estimated Trajectory of AUV’'),grid 

xlabel(’ X (in feet)’),ylabel( ‘Y (in feet)’), 

gtext (‘Estimated :Solid line’) 


gtext (’Actual :Dashed line’), 
hola off 
end * est opnp 
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(zhO,P,data, lambda, prms) 


(zh,P,2Z2s])=scen curl 
function applies Kalman filter when there are 


Punee 1 oOn 
* modified April 2nd 1994 


* this 


TS 


n=length (data) ; 

VAS, Lj arala\0) 

R=1; *Measurement error Covariance. 

al=prms (1) ;bl=prms (2) ;b2=prms (3); bl=prms (4) ; Ts=0.0225; 
FOr we — ia 


rho=deatait,2) 
alpha=datatt,1); 
headh=zh(4,t); 
RPM(t)=data(t,6); 
RUDDER (t) =data(t,7); 
a=cos (headh*pi/180) ; 
b=sin(headh*pi/180) ; 
A= [0 0 “a00 0" fo 


f 
O; 
; 


0; 


be 
B=[(0 0:0 "C7a2 050 0-0 b2-050 oO 2ere 

[Phi,del] =c2d(A,B,Ts); 
dx0=rho*cos((theta+theadh) *pi/180) ; 
dyO=rho*sin((theta+headh) *pi/180) ; 

x0=2h(1,t)+dx0; 

yO=z2h(2,e) +ayo; 

ZS (yt) —lscOy 0k: 
[v,dvx,dvy]=vp_open2 (x0,y0, lambda) ; 

h= (dvx,dvy,0, (| -dvx*dy0+dvy *adx0) «p17 1307070 0s 
S=h’*P*h+R; K=Phi*tP*h/s; 

e=Q0-v; 

zh(:,t+1)=Phi*zh(:,t)+del* [RPM(t) RUDDER(t)]’+K*e; 
P=Phi*P*Phi’ -K*s*K’ ; 


eo} 


SS 


a 


SSS eS 
SE SES hk Oa a) 
SD HS Tal Sb eS | 

OO 17] oo 
SO fF ea = 
S45) oe 


end 


3 


engys Scan cur - 


002900000 
000060 O 


2. 
° 
zhm, zsmt 


data file must contain: [bearing(relative) , range] 
estimated states are zhm=([x, y,v,thetal] 
all angles are in degrees ranges are in meters. 


0000000000000000000000000000000 
BOTT GOFTOOTFT FDO OFF FFT FOO FOTOEOHO OS 


S 

1c 

4 

* Modified February 3 1994 
% 

% 

% 

COCCCC CDEC EDEEEEBECEES 
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0000 
6600 


o\? 


! del tezgl.met 
eug; Nnold off; 
subplot (221) 
n=length (data) ; 
me—200 ; 

fer i1-=0:n/200-1 


[zh, Pklk, zs, xsmth] =scn_act2 (zh0O, Pk1k, data ((ns*i+1): (ns* (i+1) 
ie); Lambda) ; 

P—Zo ls: ) peym=cei2,:) ; 

Zam= | Zam, zh): 

zsmt=([(zsmt,xsmth] ; 

showl1l (xm, ym, lambda); title(’POOL’); 

1£ j== 

meta 

end 

Mora. ort 

lambda=0.5*lambda; 

Mmiz—length (zh); zh0=zh(:,nz); 


end 

Mieka off 

show10 (xm, ym, lambda) ; 

Mena On 

meee zaimi{l,:), zhm(2,:),'’*g')}, title(’ Estimated Trajectory 
of 

AUV’ ) 

meta off 

show10 (xm, ym, lambda) ; 

mmeoerzsmt (1,:)},zsmt(2,:),'*r’),title(’ Smoothed Trajectory 
@QmuAUV'), 

mela “off 

meta 


° 


end * est _act2 


999999999°9°9999999999999999999999999999999999999999999999999 

ooo oo O70 OC OO 0 0 COCO COCO COCOCOCTO CO OO OCCOUCOCUCUCUOCClOUCOCUCcOU CO OC COCUCUOCUCOCUCOWC OUCOCcUOCcUOCUcUOWCcUoOCUcUcOCUCcUOCCcOCcoOCUcUcOCUcCUOWULCOOUCOUCOULCOUMCOUCOUCOUC8lUL8 

Umetion [zh, Pk1k,zs,xsmth]=scen act2(zh0, Pk1k,data, lambda) 
[zh, Pk1k,zs,xsmth]=scn act2(zh0, Pk1ik, data, lambda) 


Meaified Pebruary 3rd 1994 

The estimation algorithm for navigation of AUV. 

This function works with actual data taken at Pool. 
Fixed interval smoothing is applied to the estimated 


OP off of© olf olP off olf Fh ol? 


states 
CCCEEEEEEEEEEEEEEEEEEEEEEEEEEESEEEEESEE SES ES ESS ES ESSE EECSESSS 
n=length (data) ; 
Za: , 1) =zno:; 


iS 


at=0i20225.; 
k=]: 
for t=) 


headh=zh(4,t); 
rho=data(t,Z2); theta=datale, 1); 


A=[0 0 cos(headh*pi/180) 0; 
0 0 sin(headh*pi/i8s0) 0; 

Oat OOO: 

Oreo) 10 eh) ee 

B=([0 2070 01 

[ 


Phi del] =c2d(A,B,dt); 
dx0=rho*cos ( (theta+headh) *pi/180) 
dyO0=rho*sin( (theta+headh) *pi/180) 
xO=Zih Gl? ©) eroece - 
VO=Zh (27 a0): 

ZS | = eeO yr 

[v,dvx, dvy] =VP(x0,y0, lambda) ; 
Kalman filter 

= (Ve Givi Orn. -dvx*dy0+dvy*dx0) zo Oral tsO| 
KG=Pk1k*H’ * (H* Pk1k*H’ +R) *(-1); 
Pkk= (eye (4) -KG*H) *Pk1k; 
Pkkc=[Pkkc Pkk] ; 
Pk1i1k=Phi*Pkk*Phi’' ; 
Pkikc=[Pklikc Pki1k] ; 
Zhi(:,t)=zh(:,t)+KG*(0-v); 
Zale (27) =Ziwet ae) 
2n( 2, t+) =Paie chaste 
ZNC=|(2ZNe Za nes ale 


end 


6 Smooth data points 
[xsmth] =smooth!zhc, zhic, Pkkc, Pkikc, data) ; 


end *scn_act 


CESEESESESEEEEEEEESEEESEEEEESEEESEEESEESESSESSEESESE ESSE ESS ESSSEESS 
function [xsmth] =smooth(zhc, zhic, Pkkc, Pklkc, data) ; 
function [{xsmth] =smooth(zhe, zhic, Pkike Pielke, datas, 


Modified February 11th 1994; 
This Subroutine implements the Fixed interval smoothing 
to the previously stored error covariance and state 


estimates. 

This function works in reverse time 
EEESESSEEEEEEEEEEEEELEEELLELEL ELE LL EL ELELEEE EEL EL EL ELE EEE EES 
Flip the Previously stored Matrices 

kkc=f£liplyr (Pkkee 


WS 


Pkikc=fliplr(Pkikc)  ; 

Mee=fliplr{Zhic) : 

mafe-plipilr(zne); 

Sinitial smoothed estimate is final estimate zh; 

Mere: , 1) =znc(:,1); 

siMnitial smoothed error covariance is final forward time 
6 error covariance 

Patatliplr(Pkkc({+,1:4)); 


fOr 1 =1:length(data) -5; 
headh=xs (4,1); 
A=[0 0 cos(headh*pi/180) 0; 
O 0 sin(headh*pi/180) 0; 


02-0 
OOD 

OO 4020)" > 

fentdell=ec20 (A,B ,dt): 

¢ Smoothing Filter. 
Aksfliplr(Pkke(:,4*1-3:4*1))*Phi’*... 
(fliplr(Pkikc(:,4*i-3:4*1i)))*(-1); 

pete +)1)=Zhic(:,1) + Ak* (xs(:,1)-zhc(:,i)): 
Peete laiplr(Pkkc(:,4*1-3:4*1))+... 
Ak* (Psl-fliplr(Pkikec(:,4*i-3:4*1i)))*Ak’; 


0 0 
0 0 
B= [ 


end 


6 Flip the smoothed estimates to forward time 
xSmth=fliplr (xs) ; 


end *¢ smooth 


o\° 


CESSES SCC CESSES EESEECESEESECEEESEEESEEESEESEEESSE ESSE SESS EESESES 
rain2.m 
modified December 17th 1993 
This m file applies momentum with adaptive learning rate 
algorithm to train input vectors. 

Imput vector = [x0,y0] location of the sonar tip 


ae | 
NO of? 


OP a9 of9 ofP oP olf olf oO\P ol? 


Output target vector =[v dvx dvy] potential function and 
its derivatives 
2-99-99-929299999-9-9-9-9-9-9-9-99-99-99999-999999-99-9999%9999999999999999.9 
Coooooogc ceo coc CoC coe eoeoCeocCUlcOCUlcUOClcOClcUOCCcOClcUOCcOCUlCOOC OCC OCmUCOCUlcUCcUOCClcUOCeOCUcCOC oOCUM OCUcCC OCU OCC OCUCUOCKOCCOOCCcCOCCcOCCcUOClCUOCCcOCCcOCCOCCOCUCCCcOCUCUOULUOWULcOUlULcO8lLcO 

ole 


*load nrin2.dat 
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+ lLoad@enroutz dat 


P=nrin2Z 

T=nroueEZ-: 

% initialize network 
[R, Ol =s2ze Ce) 

Si oe 

[S2,Q] =size(T) 
[W1,B1]=nwtan(S1,R) 
[W2,B2] =rands(S2,S1) 
6Training parameters 


disp freq=25 
max epoch=400; 
err doal=.002, 
momentum=.95; 
err ratio=1.04; 


display frequency 

maximum iterations 
error level 

amount of momentum 
error ratio 


of off a\P olf oP oP 


Ihea. 5 learning ratio 
Ne balers 105 increment 
lr _dec=.7 SCSCCrECMeME 


¢training parameters 
TP=[disp_freq max_epoch err _goal lr lr_inc lr dec momentum 
err Weaeoe 


io) 


% Train Network 


[W1,B1,W2,B2,epochs, TR] =trainbpx(W1,B1,’tansig’ ,W2,B2,’'purel 


layer 


¢2¢200000 
TiS -) 


poy Wy abe) 
end 
&22222922922829F2228992828999898998999999989999298929F8F89289F98299989998 
ooo o0o0 008000008 U8 U8 UO UO UO UO UO UO OC OO OO “OO OO OO OO 0 ‘oO 'o ‘oO CO ‘Oo OO OO OO ‘CO OO OO OS OO 6 OO OO 
function zhm=estneuro (data, P,zh0,prms, Wa, Wb, Ba, Bb) 
* function zhm=estneuro (data, P, zh0, prms, Wa, Wb, Ba, Bb) 
*¢ modified 11 Nov 1993 
sthis function implements neural networks for estimation. 
* Wa,Ba are weights and biases for first (sigmoid) layer 
% 
o 
° 


Wo,Bb are weights and biases for second (linear 


7° Ceoeoeeoegeooo Koo OO 
CESEEEEEEEEESEESEEESEEESES 


° 
06 CT 6 


o\o ~ 
0? 


ooo 
66% 


oP 


o$ca36 
6ES%S 


ole 


° 
x) 


o\? 


ecooeoeoeoeono ooo O88 
CECT CCCCOCCCCOCS 


Clo; Wolds ort. cite 
Subplot (2) 
n=length (data) ; 
Ts=0.0225- 
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(zh,P,zs]=scnneuro(zho,P, data, prms, Wa, Wb, Ba, Bb) ; 


aes (10 20 0 25) ) 

Meerwopn2 (25,,1 ),hold on, 

Mame Zim: )} 4 szhm(2 ej) eG. data(:,4),datal:,5) 2722"), 
title(’ Estimated Trajectory of AUV using BNN’),grid 
xlabel(’ X(in feet)’),ylabel(’ Y(in feet)’), 
gtext(’Estimated :Solid line’) , 

gtext (’Actual :Dashed line’), 

end 6 estneuro 


N 
ay 
rr ol? 
N 


modified November 11th 1993 

this function applies the estimation algorithm using 
back propagation neural networks 
potential function is replaced by neural network. 


Ph © FON $3 oP ofP of ofP of9 ol? of? Fh ol? 


Wa,Ba are weight and biases for first (sigmoid ) layer 

, Bb are weight and biases for second (linear) layer 
CCCCCC COC CCC CC CCCCCECECEEECEECEEECECEECECEECECCECECECECEECEC CESS 
=length(da Cae 
on ee, kb ) =Zh0 
=1; ¢Measurement error Covariance. 
1=prms (1) ;a2=prms (2) ;b2=prms (3) ;b1l=prms (4) ; Ts=0.0225; 
Oewe=i:n-1 


mnio-adata (ty, 2): 
alpha=data(t,1); 
headh=zh(4,t); 
RPM(t)=data(t,6); 
PUbpER (t)=data(t,7); 
a=cos (headh*pi/180) 
b=sin(headh*pi/180) 
A=(0 0 a0 0 

Cr+ On°0 

Oy alko0e 0 
9s Sas 

er O-n0 <=] 

=)0 070° 07a2 0:0 0:0 b2]; 

(Phi,del] =c2da(A,B,Ts); 

adx0=rho*cos( (alphat+headh) *pi/180) 

dyO=rho*sin( (alpha+headh) *pi/180) 

mO=Zit(i,€)+exo: 

we=Zh(2,t)4+dy0; 

a4 2 t= (KO yO]; 

* implementation of neural network. 
Aa=tansig((Wa*[x0;y0]),Ba); *’ Sigmoid neuron layer 


noo OC © 


Ge, 


Ab=purelin(Wb*Aa,Bb) ; * linear neuro layer 
V=AD (1) > dVx=AbD (2); Ovy—Ab (3, 6 network output 
h=[dvx,dvy,0, (-dvx*dy0+dvy*dx0) *pi/180,0)’; 

s=h’*P*h+1; K=Phi*P*h/s; 

e=0-Vv; 

zh(:,t+1)=Phi*zh(:,t)+del* (RPM(t) RUDDER (t)]'+K*e; 
P=Phi*P*Phi’ -K*s*K’; 


end 


end *+scnneuro 


PRDATA6 .M 
odified November 18 1993 
is program Simulates the motion of AUV in the rectangular 


OP off ol© olP off olP of? 
ct 
DS 


several obstacles are added in the environment. 
EEEEEEEEEEEEEEEEEEEEEEEEEEESEEEEEEEEEESEEEEEEEEESESEEEEESESESS 
cleamMmeigecilc 
(i= 12e ei 2 
TrS—0 02 5. 
TL=7*9; 


kmax=Tf£/Ts+1; 

x=zeros (5,kmax); 
M(:, 2 )e (2 945 20-0 
Eime. 1) — oO. 
RPM=150*ones (1, kmax) ; 
RUDDER=zeros (1, kmax) ; 
s6dynamic parameters 
al=l> 

e410) See t 

O2— 

ee 

for k=1:kmax-l1 


time (k+1) =time(k)+Ts; 
tetha=x(4,k); 
a=cos (tetha*pi/180) ; 
b=Sin(tetha*pi/180) ; 
A=(0 0a0O0O 


050 =D 2Ome 

0 0 -al O O 

QO. 07 Goer BL 

Om OF Oa 0: sales 
B= [0770-0 0 az 0-0) 0-0 ba 
E=eye (5) ; 
[phi,dell] =c2d(A,B, Ts); 
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fons del? i=c2a (Ave, fs). 
1£ time(k) >=10 
RUDDER (kK) =3; 
end 
ie tame (kK) >=35 
RUDDER (k)= 0; 
end 


Banat’ normal” ) 
ex—001* rand ° 
ey —-0701* rand; 
ev=0.01* rand: 
et=0.01*rand;: 
Sra—0.0i*rand - 
x(:,k+1)=phi*x(:,k)+del1* ([RPM(k) RUDDER(k)]‘'+del2* [ex ey 


ev 
et etd)’; 
head (k )=x(4,k); 
shdg(k) =rem m (head (k) + 0.9*k, 360); 
shdgi (k) =rem(0.9*k, 360) ; 
if shdg(k) > 180 
shdg (k) =-360+shdg(k) ; 
end 
if shdgl(k) > 180 
shdg1 (k) =-360+shdgl1 (k)+0.01*rand; 
end 
angl1 (k) =180/pi*atan2 ((3-x(2,k)),(9-x(1,k))); % obstacle 
23") 
ang2 (k) =180/pi*atan2 ((9-x(2,k)), (3-x(1,k))); % obstacle 
Bio, ) 
eee a) =rangen(x(1,k),x(2,k),shdoi(k), bi, b2); 
if (shdg(k) >=angl1 (kK) -3) &(shdg(k) <=angi1 (k) +3) ; 
dist (k) iy ate 9 ae ere x(2,k))*2)+0.01*rand; 
elseif (shdg(k) >=ang2 (k) - 2)& ee cais )<= ang2 (k) +2) ; 
Gist (k) =sqrt ((3-x(1,k))°2+(9-x(2,k))°2)+0.01*rand; 
else 
Gist (k) =r+0.01*rand; 
end 
end 
momo ik;,.)=(shdgi(k) dist(k) x(3,k) x(1,k) x(2,k) RPM(K) 
RUDDER (k) . 
x(4. kk) 5. kK) 
end 


'del datalpS.dat 


6 Save data in ascii code 
Save datalp5.dat temp /ascii 


81 


o\? 
o\° 
oe? 
oi? 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\? 
o\° 
o\? 
o\? 
o\? 
oe? 
o\° 
o\° 
ol? 
o\° 
o\° 
o\° 
o\? 
o\° 
o\? 
o\° 
ov? 
o\° 
o\° 
o\° 
ol? 
oe 
o\° 
o\° 
o\° 
o\° 
o\° 


o\° 
o\° 
o\° 
o\° 
o\° 
ol? 
o\° 
o\° 
o\° 
ol? 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
o\° 
ole 


[r,h] =rangeh (x,y, theta, Li, L2) 


[r,h] =range(x,y, theta, Ll1,L2) 
this function computes sonar range rho 


and its gradient h 


EUNCEToOM 


o\? 


02200000000000000000000000000000 
COGOCOO OOO SCO SDCO SCO DCO GCDGOODODDOODDOOS 


neta 
%EE% 


TT ol? 


my ax. 
000000 
EESEEES 


T ole 


~ oe 


° 
° 


(x,y) with heading theta 
M4 


in thevbl sy9i27 pool: 


at position 
a 
h 


OP o\© olP ol ol? olP ol? 


. 
, 
. 
, 


° 
, 


sin(theta) ; 
th2) &(theta<thl+2*pi) , 


ehays 


theta-2*pi; 
theta+2*pi; 


(theta >thi+2*pi), 


theta*pi/180 


atan2(-y, -x) 
(theta < 


aban2(-y, Liaaws: 


atan2 (L2-y, -x) 
aban2 (e2-y .o- ae: 


theta 
theta 


cos (theta) ;s 


theta 
eloull 
Enz 
wate 
th4 
while 
end 
while 
end 
C= 


(theta>= 


abe 


. 
, 


[-x*s/c°2, -1/eEnO)]); 


r=-x/c 


h 


end 


° 
, 


—~ wf 
N ~~ 
Ac rn 
LW 1 
V ~ 
1) 7) 
[we ~~ 
oD) N 
Ae < 
‘S) ” 
— ~ 
i «~ U 
— «x 
oO —— 
GQ—~N 
an 
' J 
NN WS 
C4 — 
i 
Miu i 
to Se. 
re) 
uy 
Ho 


end 
Lt 


th4) &(theta<th3), 


(L1-x)/c; 
h= [(Li-x)*s/c*2,-1/c,0]; 


(theta> 


thl+2*pi) &(theta<th4+2*pi) , 


r=-y/s; 


h 


IG 


end 
Hla 


(theta> 


[-y*c/s*2,0,-1/s] 


end strange 


end 


6 % 
open.m 
his program simulates the motion of AUV 


in a environment without borders 


smodified 25 oct 93 


° 
© 


1g 


o\° 


p 
a 
a 


OP ol? ol? ol? 


W) ole 
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cleaw,clajelc 


me—O0 0225: 

i= /*9-; 

kmax=Tf£/Ts+1; 
m=Zeros (5, kmax) ; 
See. 1) =(3.3°3 45 0 ]*; 
time (1)=0; 
RPM=150*ones (1, kmax) ; 
RUDDER=zeros (1,kmax) ; 
al=1; 

a= 00 1: 

job 

b2=1; 

BOG K=1:kmax- 1 


time (k+1) =time(k)+Ts; 
tetha=x (4,k) ; 

a=cos (tetha*pi/18s0); 

b=sin(tetha*pi/180) ; 


( 

A=(0 0 a0 0O 

GO pe<0 20 

U0 sal, 0: 0 

OF 0F 0 CO 2 

One O 0 =a) > 
Bor 0- 0:0. a2 © 0 0-0 p2i 
E=eye (5) ; 
felis delljee2d (A, B,Ts); 
[phi,del2]=c2d(A,E,Ts) ; 


memeame (kK) s=30; ; 
RUDDER (kK) =10; 
end 
fevtame (k) »> 
RUDDER (kK) = 
end 
rand(’normal’ ) 
We=Oo2Ol*rana* (i. 1 1 2 1)" 
x(:,kK+1)=phi*x(:,k)+del1* (RPM(k) 
head (k) =x(4,k); 
shdg (k) =rem(head (k) +0.9*k, 360) ; 
shdgl1 (k) =rem(0.9*k, 360) ; 
if shdg(k) > 180 
shdg (k) =-360+shdg(k) ; 
end 
me sndagi(k) >= 1380 


=33 
GO; 


shdgl1 (k) =-360+shdg1 (k)+0.01*rand; 


RUDDER (k) ] ’'+del2*wW; 


end 

ang lik) =180/pi*atane ((10-%(2,k)),,(5-x(1,k))); Ss Lee 
obstacle at (5,10) 

anig2(le = UO0y/ Di satane( (5-% 2.) jail 0-x (lek) )) 7% 2G 


obstacle at (10,5) 
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ang3 (k)=180/pi*atan2((5-x(2,-kK)) tea eee % 3rd 
obstacle at (18,5) 

ang4 (k) =180/pitatan2 ((10-x(2,k)), (2-x(Giykweee 4th 
obstacle at (2,10) 

ang5 (k) =180/pi*atan2 ((20-x(2,k)), (8-x(2) epee Sea 
obstacle at (8,20 


if (shdg(k) >=ang1(k)-2)&(s 
dist (k)=s@rt ( (5-x(1,K) 
elseif (shdg(k) >=ang2(k ‘- 2 


hdg (k) <=angl1 (k) +2) ; 

(oe) 

) 
dist(k)=ceame ((10-x (1 7k) 

) 

) 

} 


+(10-x(2,k))°~2)+0.01*rand; 
(shdg (k) <=ang2 (k) +2) ; 
“2+(5-x(2,k))*2)+0.01*ranee 
(shdg (k) <=ang3 (k) +2) ; 
“2+(5-x(2,k))*2) +0.01* vane 


& 
elseif (shdg(k) >=ang3 (k) - aa 
ees k) <=ang4 (k) +2); 


2 

dist (k)=sqrt ((18-x(1,k 
elseif (shdg(k) >=ang4(k)-2 
dist (k) =sqrt ((2-x(1,k)) 
elseif (shdg (k) >=ang5 (k )-2) 
dase (kk) —senre (Ss —x. ig) 


+(10-x(2,k))*72)+0.01*rand; 
eerie isa 
+(20-x(2,k))°2) +0.01*rand; 
else 
dist (k)=500+0.01*rand; %* add very large number 
end 
temp (k,:)=(shdg1(k) dist(k) x(3,k) x(1,k) x(2,k) RPM(k) 
RUDDER (ke) 2564. ke) (Selo 
end 
del datalpi. dat 
% save data in ascii format ina data file 
save datfebl.dat temp /ascii 


84 


LIST OF REFERENCES 


Friend, John R., Design of a Navigator for a Testbed 
Autonomous Underwater Vehicle, Master’s Thesis, Naval 
Postgraduate School, Monterey, California, December 1989. 
Zyda, Michael, and others, " Three Dimensional Visualization 
of Mission Planning and Control for the NPS Autonomous 
Underwater Vehicle," IEEE Journal of Oceanic Engineering, 
Pete 154 NO 3 pe 1217, July 1990". 

Burl, J, B., "Derivation of Kalman Filter Equations, ” 
notes for the EC 3310 (Linear Optimal Estimation and 
Control), Naval Postgraduate School, 1990, (unpublished). 
Meditch, J. S., On Optimal Linear Smoothing Theory, 
Technical Report 67-105, Information Processing and Control 
Systems Laboratory, Evanston, Illinois, March 1967. 

Ljung, L., System Identification Theory for the User, 
Prentice Hall Information and System Sciences series, 1987. 
Percin, E., Sonar Localization of an Autonomous Underwater 
Vehicle, Master’s Thesis, Naval Postgraduate School, 
Monterey, California, 1993. 

Cristi, R., "Sensor Based Navigation of an Autonomous 
Underwater Vehicle," paper presented at the International 
Symposium on Unmanned, Untethered Submersible Technology, 
Durham, New Hampshire, August, 1993. 

Beale, M. and Demuth, H., Neural Networks Toolbox, The 


Mathworks Inc, 1992. 


85 


Or 


i. 


INITIAL DISTRIBUTION LIST 


No. 


Defense Technical Infemnatkaon Center 
Cameron Station 
Alexandria, VA 22304-6145 


Library Code 52 
Naval Postgraduate School 
Monterey, CA 93943-5101 


Chairman, Code EC 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5121 


Professor R. Cristi, Code EC/Cx 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5121 


Professor H. Titus, Code EC/Ts 

Department of Electrical and Computer Engineering 
Naval Postgraduate School 

Monterey, CA 93943-5121 


Deniz Kuvvetleri Komutanligi 
Personel Daire Baskanligi 
Bakanliklar, Ankara, TURKEY 


Deniz Harp Okulu Komutanligi 
Tuzla, Istanbul, TURKEY 


Golcuk Tersanesi Komutanligi 
Golcuk, Kocaeli, TURKEY 


Taskizak Tersanesi Komutanligi 
Taskizak, Istanbul, TURKEY 


Alp Kayirhan 
Lise Cad 92/4 
Denizli, TURKEY 


Ziya Kayirhan 


Lise Cad 92/4 
Denizli, TURKEY 


86 


Copies 





a i, a 


Or 


' ae : ‘ _ 
eet! - . er 7 


- 





DUDLEY KNOX LIBRARY 
NAVAL | KOVATE SCHOOL 
MONTEREY CA 93943-5101 














rere OT 
sein iat iat dy Hea ro te a | DUDLEY KNOX LIBRARY te Seth act ee cree 


‘ , 
Bai 
CRSP ORNL Ppa feiss ee cagtacict ee f ie bhth «al mT rin MUTATE TT 
Aste fyte hs Sag Ei Set ah H £4 Tact ugh egtny f 4 i ‘ Wy a | ! | | HE a a a 
cat i B ries sage Ny he rH slat f t J | Waal || | 1 ALA aba II 
pt 4 wi 4,9) LA, b 4 rf abe’ ad 44 git: }H a | | 
it deep can rab apes pera matue task pest Gite stay SME att | WANA | 
Sebel petite etek ete Life pees AME aig to il MMi PL SUE 
trp dales Dota & bes 5 } eget yh : Mei, 
i beat ti: : ¥ ie 2 abt ante ie ia Mitte 4 < 8 
Bah hn Lae of, eaes it, Ue ts bal . > kia rath : j > fehetce 
Ms 
- 


Drfogeyse 
a EY 










te wei jo. s 
| Picea ded SEE 


3 2768 00309493 9 








ey 
ais 





























ake ‘D> Ae ‘ Tatas do nt wor t. ‘ 
vee . i a ‘ yal es | Pheer aie ‘ ae vas * 
re Has ihe wth iantiets (iia? (eet? trae! pO) | TE oes ‘, epeee 
% etal ls, 7k Petste ils eo SEO “a te fede meg oe cer eer AL Wisse _ 
Arelal o pemesice weg! Pa SUR Ie Pubs, 6 Sccreiee aimee se er 
eee | car a7) 698 Re ee ; > 
PRT Pt ae : H Wee \' 
r 45 : ie aan 4 ues yqucea tee 'j 
, hee Lh a fastens o 1 dere 
ty * (is { 
Ve tebe 4108 t 
e ri aie fer} 
vee 
tone 
sareefye! 
i 
64 ¥ ee RF ogee 
TFs, i 4 Pee: , Rel 
: X ¢ . is f , ’ ’ 
ed ate? } U ‘ t z ‘i ; Wa. 
ay 5k 


e Pye, 































ceheine 
Aa chet ld Wig i mone 
13%) 0.4 ae bs he = ‘ 
gla po dg sige. : Neh k ee i PA bh eek hehe slob ' ! -} fa Le sie 
Pi Pad, 5F 3 fay B . » 8 ant { ke ada 
Api bt ‘ ihe j eatetuss 
ty hb dys fed 
wat Cartes} TY ae sy ” a 
: ‘ 
ay ; "stetae ¢ ries 2% ere 
: rhea hager 8 i. Ps a 
heel ab Tat pe * 
H tal off org saat atta : ' 
j ¥ . ‘ ‘ taayt ‘ 
ib ans cag Spgey bya! H+ SctB. Vib el wleteeciesine . 
; cir er decnt af ae :! Poe ad Fibed Soh every = tte ' > ; ’ 
. 
i ated és, § boa x !,, Pyles? Glin Bass a) | t tt ; at 
i hae ph A ver Savane Pag eos Teor. ent : i , ' 
Fy ALU eae ed ole? ’ ofy tht bt, beeed pyar 6 1 We eer mg 
Pp 2 3 * Behe N tee hee Sat Cee ature ae eo tees tane : , 
hee pe 3+ ie he <8 ost %on peek oe eno ec Atos & 
rane Hue ey Aer ee Lied foot! rs 
oe eur. ey. 5 3° os 
i oe dat e lee a ae ahi 
1 a eee eR Ar Bat eica 
at me Ao ene a eu ire he 
ety, ae Pad Aa cet al te ar] Res 
ald 






aM pecs 
Fahy s, 
BAe as BPSS. 
Nee Pad ras i 
JOM A 4 
hh 


ae Ub 
















WG fo 










+A a bet , 2 ame 
HOU OM ie ie 
' = t. 
ait 
# trparte r) te 
Pisenhte Pet & ape Pian is vine 





5] Bea Raat ea rd 


on ea? Ua 






TY Fe 
vial 


=e 
“"e 


Poe Palo Fe 
ory 
. 
af stye 










tea 
eats OFF 8s ag 

Yo tas 
“ECE Wt dae al 












1 
HP Nee 
+ Fee SPD tee 
Veet Mete tig? eee 
* aozee gad 








ve ee 
Pt get 


‘i ee Tn i 
Merepert ate ov Oe aL HL J ’ 
a a) * 





Sr 










" 
at FPL Bye 


Polit doe 
ene e athe 









ow we a 

Se ue =, 

Pew 
ee ed 
am 407 * 





t try nay 
LLY ae oe 
‘ 





setheber 


‘7 t 
PNP a Peay 

















Ayrbebarege 
‘ ‘ 








37 





~. 


+ ene st we) YE were owe 









wo ae 
Fok 
= 








yt 





eo 


90° btm 


ed 
ee ae 






‘ 
beh vetG Rae 
181 ae dhe et ohane 8% Rbers 





a ae Bee 
re a 















— 





= wae a 
em ar Me 

















































































V mat hy per ggeed Ygpwe Mere uv 
Hi b] hee) ee 2 9 UA Fab 
3: i Sov plete Rr Cher ae) 
Hi a A be TE: : eae 
i aent ony ten t ef a 
: olay ’ a0 : a 
i Py wi ' eo» a 
a Fass ghtirogfertving yt Pett ogee ceed, ge to yk ee : 
{! Che FRAC MOO a red bee aks AO SR cant a, arate yg nr 
i Foe Ps4t Fusgerlipe CUS, ve yet gehen gal dove tede 1 
% iF Toph ebes yr ys a gttee abet done oy aon ow bY oe ait ae 
“es ys ate O5hy Tap Oy at 1 a aha pte aa. 4 ‘ . 
: # F3 1h. > abe *. 
s 4 } sbeps 
t were ater 9g, 4 © 
? cb attate gs ; t 
4 AA 
hi Tt beh Phas Were n 5 ’ 
$ Un PUT set eae de megs ““se veel t ese 
5f Aste tt nas | er aN of ule hee 
, 5 wee a 4 “. fod4 oe 
Hi ‘ Payee? Bar ieenorrsp Ue alenee mete te tt ut 
5 Sox ety OES OG meh s yeenr eee be" % 3° 
}, TUT aria Sm et ta ed 3 j teed, 
3 ot zOeus et el cee Hay eel ee $7 Ge 
cad weer ayenp ayes A wipe mates «8, he, cry hays . went ap 
4 fp, See ght yet st ie = ' PVF tet eh cy ‘ ‘ eof 
seas be Uetishserics + t x ithe Ri Pteed Bh Ge bah, vy i DO s it # «pn ey 
a | UT PO pgs gre 1s aM iad bbe) teed : be FF reek tm Pasery aves oot ay 
. 2 vs rhea Pe ben | teehee? » re Braye H 
ant ater y. Wee bee Phy ee +cblase® ohedks Pure eebare Qsteny, 
iz ‘ een td Lae 4 . ‘s a ed ed ee © 





Jor Aya pe,’ 4 
enn eel f 
Peat otatgt 


em we 





2 iS 
; 


Vetefbtee epee oe | 
ete 1 
¥ 


¢ Ro neae 





w~> 







gaa 

By i eseconssste maha | 
0 1.) topes be feet, a 

PSS SPS elas Pe heard oF 
a’ phe FF ahd be 2 

4 






on 


ae “yJW aby Ppe tgs, 


r 
Res ’ i 
wipide se ay?ss ala 























































































































































































































































































































































Wats A HM ; ‘ . ' ul ‘ ? ; : Pasay Y Leite og 
aids Rb VE atin aap a Fayy« Hin : hes Fars AY SUH AEM Hey " SNE , 
} 7 K H Tel ul ty gey deg eahersta y Ce Fee 8 ed | oon i te he Bes i beh 9 COL 
PAYS S glad tenets pte fae e Pee Pi Se LY LY Pre ah on ty T8084 bse cis) ppareh SINT 
rane dene?” es to ophyg en of} ‘ 5 : sere es " fa gitewy, ey i (aN, os s Pe bin “ 
peur oy ’ . “ Pa roo ef . ' 

Ropes theve eary pone ec Apr Pe eaten te eatey sila Star ate ' wea core ais) ow Ne lea eee ty GRY 
paseo susigntg pee WPS NR Ea aad Me aqes kites bags 
ane - ‘ 2 6B ae oF on Ut > ew 4. 4p wae ‘ ‘ ; 

ies ar one “Pea et Gulab : a ats t ; idivee na meee een fa ate 
att oe Cored Mt ae OP san ret® Sheep pod * ‘ U deus ’ aay ee 1 
= | ‘ te eles ose 
cate. Ms Sie Ip hay co aoe aca beae te tes yi ‘ _ mehr ae aaa 
? ‘ : ‘ 
Debye Lead ra « Fanaa oe peta. 
As) AL Sd at ’ wize oh in ebtse Hae Peep fe gt hy ‘ x ’ ue - ; 
in Li t + WerseAPbsete gs Ws eyeiebas bonote dob 1 ue Aaa eres Es iat LT) 
Bey 43 ayy? Lior ore SH te ay th = 
¥ a > 5 Faia? PPPS sop a tay, brhte want ‘ ' a0 . ee F 
rt y Seip : . at ODEs Sowa ty logy aan 7 Atte seve ate 2 ; ae 
: r es 2 wee i ’ - 3 4 e4 4 . . ] aie ' 
6 OOP oF gear, ott Sas ta 4: 3 ; : ’ 1 geese a : HEAT APP eae oe Mesh asuir ales dye) ae CEE 
IAP t dere 4 Fee * : : Porter) A seu Set wastes ‘ SLE UE Boe Clb, 
> 4 ¢ vere a *Oo* mata 
1 0 8s Leh Ae re he tae Gout dle Pay er sari 
PS at bab ed A pu Bt te ly take 08 w i, ' 
ess’ TS ‘ t ote Four a tg pe ttin ys Ce an rahe 
- * k aL Uae of, g* e 
postage ey ve * : ; : 4° we 
oat baiplers 7 Dera atte, une PU eWeb ete the Petit a0 ve Set 
PASR TE Bewhg spear f DOP ntyes peor ogg ae he teh US s cng? 19 cee t sbed Fe ’ ' . 
AD Ld Ke bye Leelee elas ober atyy Jae 13 Serer Lis Ueeriee ; ’ tei 
‘ Cede tes Sit fk yiattes te ee ts PRM ee yt Par Lata ree ta ete 
i thetic Ps ches. Shin Sha) ipeedune de® 3 Seth a Bh aes is in Uo aaa Ve re, Aid a 
Vey, Wa urine ApOge fem shi yr verueti an red g ay’ NE aT eONGS OF tae rer hre gly 1 ’ 
irpte) LEY cia, atin, 3 pS Oe Utree be we to tpacye : 
Od te 2h Tpger & be has ic EON A A O BRe cel aap 8 aac Sea drone eho 
ety Fy as reer) WAN 13 bate Wi ia oases Be oeH tepens va, Rare sons 
ob ' vente S785 tbe, ee Oe ong au, PH OE | 5 ‘ i as =n 
heared] : pes a0 Odea Pe it t ene ie : 
. oP abctasate? Dees yy ACC , i] ao cay nhs 
] US Peep Ht 4 © 4b toc agn Pr ererere i " Pir 
botante gl Eas Mn tad od ALE & 0. wy eta 
124705 of OT Rse tee A gb stmt Spe pei ys , =a { 
Ab efa Ue ote Seu pr Ore eruee® fe ba, ma tee é Petiteg 5 Uy ae tre 1 
bg ais TTR why ' OREO a RTD he oy Hist +4 eran apeebe . Pa ae bhes r ata aan 
4 ren ee et begs P ara e 4n Pat ye i J Ure Pe ee . 4 
' sapintets | ed Ke t ;" Vbedease Ae irrety les ee oor a rf mors Ui eee eal a % ee im we 
np re ee veya BG bre. hen Tn Pale | ste wae hep badat Gf Han rs toe 
40+ Dbee Pargeten Tatar ole Ie AS te eee ee ed pease pe dy vor 8 ee ‘ 
eR pBEpete te ‘ vi pSbags ted tar Le Wik bee 29's ose bes @ hey vac pk , one Ae 
Paras ae 5 water ae ig gegcal tay sett wites ei gee at, eet tetes ere hes at ea Pag ace s4 
Tire ge FOOT UBIDARS hep yt eyVlecerges omer re LPT ' At hace OP Or GES Asp egee® re geet g ray rice: 
thas ; ', rol gf te Gs ay i vet te bee 8 meee. kel pH th DALY ATT Td WY TCT a Wire! 
A f ‘ 102%, OPO Te 8 wt Bw, - COR CPL: afin Teta ttuae nie aaa ‘ 
4 Hpi ohn. fee a ote gy pcb en Witter ee ar non u 4 
ry A ‘ i 4 viey bat feteare & We och wet 6,, of elas reer ei? 
Xe hk ed or2 g $ Dette Rie ’ a et ee are See aro aie ‘ 
4 $i ee eo t { j s + \ b t ' t t. aye-4 1 wie Fi 
Si 774 He 4 may + Yoge etn ol ian os : 
ah, ® shee ore cm Eats dele a be ae wa 
oF aa Lt Hate’ ys H ¥ ral to ote ; phe % Be a! "e = Ae : a ate ot , 
ib pe Uh 3 ; : aia . FELO Vowel reap, + # th ' =r : 
eRetarwe bay icy M f : ? ¢ Pree ly Go arty SA Err eee Cree tbe SSS Weiner sero bay, 
ate swede t s Pe Peale ny yr ah. tag oc HY g ANTS crete A Piee Uren ic eee ea 
Tyrer t oe t ‘ fas gaee lew 
> S14 2 of ' 49d Somes a 
46d GAL BRB lat ae aac ta ere seth : 
pone Kvle sf kD ‘ Fit ar Seer I < 
: c ’ J t ‘ i oe | al) r ong Paes . 
4 Syte ef : » weipet 8 Pel ey ° rf = 
Fotatytrige pgs tar a i . P - Ja 4 1Sépe tee ued ert at sete suse iG 
ats Fy hihi es 4 : i "oF Mlate ong te crisis ae 
s feet a bby oe pe Tw. ong rae ! 
$ Ob aenoatt ¢ Ce Ot ere CUI 
‘ “4 ‘WP . ; ae Py Z 
, peo wht gered gig ; > = baee 
bt rane weep Pa Ck aes he Pooh al ' poe 
>i + pte ‘ “4, t “ ory ‘ 3 . 
Se Gre a DPT P eae ea Send USEPA Nel sticg tt gt oe aE 
TF ab Vote viares * Wi, gtd RAlAst oe sy ' eet Yh Gib geen ‘ 
So. weg eat Gtae Yt " Bat Dane? tphob oc dae oe ] hogs tee ont OS 4 
SES Teabys : Ot PAV G sat Gefed eto ie! teay Re ees ts © al tte 
SURE ae ie fone o eine, te tet O75 un vgs tact 
‘ tf. 44 
Rod Pat i } GT) Val 6 Babine 4 1 ae Yo ‘! x Mey 
PF Mises of POTS. 6 Valle He bt aes, Pee ele @ agels fy ebe 
a arate yes ' PEL re) ee hd eee "4 ome 8 +e aa 
; : We a Fe 
Stan Te ot Oe eee Lee tere ; pers 
tone gey 1 eo LeU ALU ecm seca gals , 
rari yee ' Hoenn Neaeeee Mega CSL Siig fat 
* RVEG 49} sp eras 7 Lae SeaClear ai rho an ‘ 7 ’ 
+f Sr ute* ou i 5 ‘ eye Pate i as z a at',* i? 
Wher it G . Mas rar conti ‘ - rou Pa RoC Hal Te mart eerie Rr ioe ene 
1°17 phase x iy etarS ff i") ' POCMLIETE  Wic vat Net ORS 1 seis wir yleerant : 
4 ue hd fi tp.  @s peepee ’ ke 
Eee ede te ue say? k aby - : 





[ei ee clon a td 2 


ye 






Woeagp 
ther sel, 
Vyas 


erfctra 


Tea Pj ales 


¥v 
PADS Be yew 
Day 
U BY 





